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Benchmark  calculations  for  higher-order  parabolic  equations 

Michael  D.  Collins3’ 
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(Received  14  October  1988;  accepted  for  publication  6  October  1989) 

Benchmark  solutions  generated  with  parabolic  equation  (PE)  models  are  presented  for  range- 
dependent  underwater  acoustic  propagation  problems  involving  both  penetrable  and  perfectly 
reflecting  ocean  bottoms.  The  solution  of  the  w  ide-angle  PE  of  Claerbout  [J.  F.  Claerbout, 
Fundamentals  oj  Geophysical  Data  Processing  (McGraw-Hill,  New  York,  1976),  pp.  206-207) 
agrees  with  the  outgoing  coupled-mode  solution  for  the  problems  involving  penetrable 
bottoms.  The  solution  of  the  higher-order  PE  of  Bamberger  et  al.  \  Bamberger  et  at..  “Higher 
Order  Paraxial  Wave  Equation  Approximations  in  Heterogeneous  Media,"  SIAM  J.  Appl. 

Math.  48.  129-154  ( 1988)  ],  which  is  a  generalization  of  Claerbout’s  PE.  agrees  with  the 
outgoing  coupled-mode  solution  for  problems  involving  large  variations  in  sound  speed  and 
propagation  nearly  orthogonal  to  the  preferred  direction.  The  computer  code  FEPE  was  used 
to  generate  the  benchmark  solutions  and  was  found  to  run  several  times  faster  than  the  IFDPE 
computer  code  due  to  a  tridiagonal  system  solver  in  FEPE  that  is  optimized  for  range- 
dependent  problems. 

PACS  numbers:  43.30.  Bp 


INTRODUCTION 

The  accuracy  of  the  parabolic  equation1  (PE)  method 
in  underwater  acoustic  modeling  has  been  assessed  with  nu¬ 
merous  range-independent  benchmark  problems. :  The 
wide-angle  PE'  has  performed  well  in  most  of  these  tests. 
Several  range-dependent  benchmark  problems  were  recent¬ 
ly  posed  and  preliminary  results  were  presented.4  Some  of 
these  problems  involve  perfectly  reflecting  ocean  bottoms 
and  thus  provide  an  extreme  test  of  the  ability  of  a  propaga¬ 
tion  model  to  handle  wide-angle  propagation. 

Since  the  standard  wide-angle  PE  cannot  handle  propa¬ 
gation  angles  much  larger  than  40  deg,  a  higher-order  PE 
model5'6  based  on  a  Pade  series  has  been  developed  to  handle 
these  problems.  The  higher-order  PE  accurately  handles 
propagation  nearly  orthogonal  to  the  preferred  direction 
and  produces  solutions  essentially  identical  to  the  outgoing 
coupled-mode  solution.7  The  computer  code  FEPE*  is  used 
to  generate  the  benchmark  solutions.  An  efficient  tridia¬ 
gonal  system  solver  (not  based  on  the  standard  Gaussian 
elimination  scheme)  in  FEPE  is  discussed,  and  FEPE  is 
found  to  run  several  times  faster  than  the  IFDPE1*  code  for 
one  of  the  benchmark  problems.  Solutions  generated  with 
the  coupled  normal-mode  model  COUPLE'”  are  discussed. 


I.  THE  HIGHER-ORDER  PE  MODEL 

Solutions  generated  with  PE  models  approximate  the 
solution  of  the  outgoing  wave  equation 


~  =  ikd  1  +  *Q, 

dr 


=  k„  2(a 


.k>+iL-±3e.±\ 

°  V  n  dz  dz) 


(1) 


(2) 


dr  p 

We  refer  to  Eq.  ( 1 ),  which  can  be  solved  in  terms  of  outgoing 
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coupled  modes,  as  PE  ,  .  The  reference  sound  speed  is  c„ 
where  k{l  =  co/ct),  r  is  the  range  from  a  point  source,  z  is  the 
depth  below  the  ocean  surface,  k  is  the  complex  wavenum¬ 
ber,/?  is  the  density,  and  to  is  the  circular  frequency.  The  PE 
field  Q(r,z)  satisfies  homogeneous  boundary  conditions  at 
the  top  and  bottom  boundaries  of  the  waveguide  and  the 
initial  condition 


(?(0,z)  =  v2tr/ 


I 


(t>l(z{))d>l{z) 


(3) 


where  z„  is  the  source  depth.  The  normal  modes  d,  and 
eigenvalues  satisfy 


F^_±dp_ML  + 

dz~  p  dz  dz 


k'A- 


(4) 


In  practice,  either  the  Gaussian  PE  starter1  or  Greene's 
wide-angle  PE  starter 1 1  is  often  used  to  approximate  Q{  0,z ) . 
In  range-independent  environments,  the  complex  pressure  P 


is  related  to  Q  by  Q~>[rP  for  /c,,r>  1.  In  range-dependent 
environments,  d  /dr  and  x  do  not  commute,  reflections  can 
be  generated,  and  a  term  involving  dp/ dr  appears  in  the  wave 


equation.  Thus  Q~-frP  in  range-dependent  environments 
only  if  the  range  dependence  is  weak. 

Bamberger  et  al.  used  a  Pade  series  to  approximate  the 
square  root  in  Eq.  ( 1 )  and  derive  the  following  higher-order 
PE: 


dU, 


dr 


1  +  b,.„x 


U„, 


a,.„  =  [2/(2  n  +  1 )  ]  sin2  [yV/(  2n  +  1)), 
bjn  =  cor  [Jit/ {In  +1)], 


(5) 

(6) 
(7) 


where  Q=tQ„  =  U„  exp  (/£„/■).  Equation  (5),  which  we  re¬ 
fer  to  as  PE„ ,  has  been  applied  to  underwater  acoustics  and 
solved  with  the  method  of  alternating  directions.6  This  ap- 
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proach  involves  n  steps  with  the y'th  step  requiring  the  solu¬ 
tion  of  the  equation 

( 1  +  b,„x)  =  ik()a  nxU„.  (8) 

or 

Since  Eq.  (8)  is  of  the  same  form  as  Claerbout's  equation  (or 
PE,),  for  which  simple  and  effective  numerical  solutions 
have  been  derived, 11  “ 1 '  the  alternating  directions  solution  is 
easy  to  implement  into  an  existing  PE  computer  code. 

The  computer  code  FEPE  solves  PE„  using  finite  ele¬ 
ments  for  depth  discretization  and  Crank-Nicolson  integra¬ 
tion  in  range  as  described  in  Ref.  6.  The  tridiagonal  system 
solver  in  FEPE  has  been  designed  to  minimize  computation 
time.  The  code  uses  an  elimination  scheme  that  involves 
sweeping  downward  to  the  row  corresponding  to  the  ocean 
bottom  to  eliminate  entries  below  the  main  diagonal  and 
sweeping  upward  to  the  ocean  bottom  to  eliminate  entries 
above  the  main  diagonal  followed  by  back  substitution 
sweeping  up  and  down  from  the  ocean  bottom.  In  contrast, 
Gaussian  elimination  involves  sweeping  downward  to  elimi¬ 
nate  all  entries  below  the  main  diagonal  followed  by  back 
substitution  sweeping  upward. 

For  problems  involving  range-dependent  ocean  depth, 
the  new  scheme  is  more  efficient  than  Gaussian  elimination. 
In  the  decomposition  into  upper  and  lower  triangular  matri¬ 
ces  of  Gaussian  elimination,  it  is  necessary  to  repeat  sweep¬ 
ing  downward  from  the  ocean  bottom  as  the  ocean  depth 
varies.  With  the  new  scheme,  it  is  necessary  to  repeat  sweep¬ 
ing  only  for  a  few  rows  near  the  ocean  bottom.  Since  multi¬ 
plication  is  faster  than  division  on  computers,  the  tridia¬ 
gonal  system  solver  has  also  been  improved  by  replacing 
divisors  with  factors.  The  code  FF.MODE*  determines  the 
eigenvalues  using  the  finite-element  matrices  and  constructs 
<?(0,z). 

II.  BENCHMARK  PROBLEMS  AND  RESULTS 

Problem  1  consists  of  three  parts  each  involving  a 
wedge-shaped  ocean  in  which  the  sound  speed  is  1500  m/s, 
the  ocean  depth  decreases  from  200  m  to  zero  over  the  first  4 
km  from  the  source,  and  the  surface  is  pressure  release.  A  25- 
Hz.  source  is  placed  100  m  below  the  surface.  For  part  A,  a 
line  source  is  used  with  a  pressure  release  ocean  bottom  in 
plane  geometry.  For  parts  B  and  C,  a  point  source  is  used  in 
cylindrical  geometry  with  sound  speed  1 700  m/s  and  density 
1.5  g/cm'  in  the  half-space  sediment.  The  sediment  is  loss¬ 
less  for  part  B.  The  sediment  attenuation  is0.5dB/T  for  part 
C 

Since  energy  is  reflected  back  toward  the  source  by  a 
pressure  release  ocean  bottom,  PE  ,  cannot  provide  the  full- 
wave  solution  for  part  A.  However,  we  apply  PF„  to  this 
problem  to  show  that  it  accurately  handles  the  outgoing  so¬ 
lution,  which  involves  propagation  angles  up  to  nearly  90 
deg  near  mode  cutoff.  PE„  is  solved  over  a  sequence  of  stair 
steps  that  approximate  the  wedge  geometry.  For  this  prob¬ 
lem,  it  is  necessary  to  remove  reflected  energy  by  mollify¬ 
ing14  Q  at  the  beginning  of  each  stair  step  as  follows: 

Q(z)^{  Q(z’)  '^d>J(z')d>l(z)dz\  (9) 

J  i  i 


FIG.  I.  The  wedge  with  pressure  release  bottom  in  plane  geometry.  Trans¬ 
mission  loss  at  depth  30  m.  The  solid  curve  is  the  PE,  result.  The  dashed 
curve  is  the  PE,  result,  which  has  large  phase  errors  just  before  the  mode 
cutoff  ranges  near  1  kni  and  2.2  km. 


where  the  sum  is  over  the  ,V  propagating  modes  at  the  range 
of  the  stair  step. 

Transmission  loss  at  z  ~  30  m  generated  with  PE,  and 
PE,  appears  in  Fig.  1.  The  PE,  result  agrees  well  with  the 
PE,  result  of  Ref.  15.  The  PE,  result  exhibits  the  largest 
phase  errors  just  before  mode  cutoff  at  r  =  1km  and  2.2  km 
(three  modes  are  excited  by  the  source).  Data  for  this,  as 
well  as  the  problems  that  follow,  appear  in  Table  I,  in  which 
CPU,,  is  the  run  time  required  by  FEPE  to  solve  PE„  on  a 
Digital  VAX-8650  computer,  A z  and  A r  are  the  depth  and 
range  increments,  zXf  is  the  maximum  depth  of  the  computa¬ 
tional  domain,  and  N0  is  the  number  of  modes  used  to  com¬ 
pute  Q(0,z). 

For  parts  B  and  C,  Greene's  wide-angle  PE  starter  is 
used,  and  the  attenuation  increases  artificially  in  the  lower 
portion  of  the  sediment  to  prevent  reflections  from  the  artifi¬ 
cial  pressure  release  boundary  at  z  —  zM .  Transmission  loss 
at  z  =  30  m  and  150  m  generated  with  PE,  and  PE,  appears 
in  Figs.  2  and  3  for  parts  B  and  C.  The  PE,  results  agree  with 
the  PE,  results  of  Ref.  15  obtained  using  the  IFDPE  code, 
and  the  PE:  results  agree  fairly  well  with  the  PE ,  results  of 
Ref.  1 5.  This  suggests  that  Greene's  wide-angle  PE  starter  is 
accurate  for  larger  angles  than  PE,.  Using  the  input  param¬ 
eters  used  in  Ref.  15.  FEPE  runs  several  times  faster  than 
IFDPE  for  this  problem  due  to  the  efficient  tridiagonal  solv¬ 
er  in  FEPE. 

The  coupled-mode  code  COUPLE  was  also  used  to 
study  this  problem.  Benchmark  solutions  generated  with 
this  model  are  not  presented  here,  however,  because  this  is 
done  in  Ref.  15.  However,  it  is  perhaps  worth  mentioning 


TABLE  I.  Data  for  the  benchmark  calculations.  Many  of  the  input  param¬ 
eters  are  identical  to  those  used  in  Ref.  1 5. 


Case 

N„ 

r>, 

A  r 

A  z 

(n,CPU„ ) 

( n.CPU„ ) 

1A 

3 

1500  m/s 

5  m 

I  m 

(3,15s) 

(1.8s) 

IB 

1500  m/s 

5  m 

1  m 

4  km 

(2,2  min) 

( 1,1  min) 

1C 

1500  m/s 

5  m 

\  m 

2  km 

(2,2  min) 

( 1,1  min) 

2A 

10 

1700  m/s 

1  m 

1  m 

t  km 

(2.8  min) 

( 1.4  min) 

2B 

17 

2500  m/s 

{  m 

i m 

j  km 

(5.1i  h) 

( 1,15  min) 
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FIG.  2.  The  wedge  with  lossless  sediment  in  cylindrical  geometry.  Trans¬ 
mission  loss  at  depth  (a)  30  m  and  (b)  150  m.  The  solid  curve  is  the  PE, 
result.  The  dashed  curve  is  the  PE,  result. 


FIG.  3.  The  wedge  with  lossy  sediment  in  cylindrical  geometry.  Transmis¬ 
sion  loss  at  depth  (a)  30  m  and  (b)  150m.  The  solid  curve  is  the  PE,  result. 
The  dashed  curve  is  the  PE,  result. 
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that  the  run  times  required  to  produce  data  at  many  receiver 
depths  (to  produce  contour  plots)  with  COUPLE  were 
more  than  ten  times  larger  than  the  run  times  required  to 
produce  data  at  one  receiver  depth  (to  produce  a  transmis¬ 
sion  loss  curve) .  This  example  illustrates  that  methods  based 
on  spectral  decomposition  can  be  inefficient  if  the  solution  is 
desired  at  many  points  in  the  domain.  Normal-mode  models 
are  usually  used  for  range-independent  propagation  prob¬ 
lems  when  relatively  few  receivers  are  involved.  When  the 
solution  is  desired  over  the  entire  domain  ( this  is  the  case  for 
matched-field  signal  processing),  however,  it  is  possible  that 
PE  models  are  more  efficient. 

Problem  2  consists  of  a  parallel  waveguide  in  cylindrical 
geometry  with  a  pressure  release  surface,  a  rigid  bottom,  and 
the  sound  speed 

c(r,z)  =  (1500  m/s)/y;  1  +  aE  +  PE 2  +  yE~'  +  SE T, 

(10) 


a=  —  (27r/i|///)cos(?rz//f),  (11) 

/?=  (nht/H)2  —  ( 4irh2/H)cos(2nz/H )  ,  (12) 

y  =  (4irA|A2///3)cos(ff'z///),  (13) 

8=(2nhJH)\  (14) 

E=  exp(  —  rrr/H),  (15) 


where  A,///  =  0.032  and  h2/H  =  0.016.  Since  the  range  de¬ 
pendence  of  the  sound  speed  becomes  more  gradual  with 
range,  we  update  the  sound-speed  profile  every  range  step 
for  r<  1  km  and  every  tenth  range  step  for  r>  1  km.  Two 
cases  were  originally  posed  for  this  problem.  However,  we 
consider  only  the  25-Hz  case  with  z0  =  250  m  and  H  —  500 
m.  Following  Ref.  1 5,  we  divide  this  problem  into  part  A,  for 
which  the  first  ten  modes  are  excited,  and  part  B,  for  which 
all  17  modes  are  excited. 

Transmission  loss  at  z  =  250  m  generated  with  PE,  and 
PEZ  appears  in  Fig.  4  for  part  A.  The  PE,  result  agrees  with 
the  PE,  result  of  Ref.  15,  and  the  PEZ  result  agrees  with  the 
PE„  result  of  Ref.  15.  Transmission  loss  at  z  =  250  m  ap¬ 
pears  in  Fig.  5  for  part  B.  The  PE,  result  agrees  with  the  PE, 
result  of  Ref.  1 5,  and  the  PE5  result  agrees  well  with  the  PE„ 
result  of  Ref.  1 5.  A  large  value  was  used  for  c0  for  part  B  due 
to  the  large  phase  velocities  of  the  higher  modes. 


30 


—  40- 
£ 


i  ■  i - 1 - 1  ■  •! 

0  12  3  4 


Range  (km) 


FIG  4.  The  parallel  waveguide  in  cylindrical  geometry  with  range-depen- 
dent  profile  and  ten  modes  excited.  Transmission  loss  at  depth  250  m.  The  — 
solid  curve  is  the  PE,  result.  The  dashed  curve  is  the  PE,  result.  V 
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MG.  5.  The  parallel  waveguide  in  cvlindrical  geometry  with  range-depen¬ 
dent  profile  and  17  modes  excited.  Transmission  loss  at  depth  250  m  The 
solid  curve  is  the  PM  result.  The  dashed  curve  is  the  PL,  result. 


III.  CONCLUSION 

The  higher-order  PE  of  Bamberger  el  at.  produces  re¬ 
sults  that  agree  well  with  outgoing  coupled-mode  results, 
even  for  propagation  nearly  orthogonal  to  the  preferred  di¬ 
rection.  Greene's  wide-angle  PE  starter  appears  to  be  accu¬ 
rate  for  a  larger  aperture  than  Claerbout's  PE.  Due  to  an 
improved  algorithm  for  solving  tridiagonal  systems,  the 
FEPE  model  is  several  times  faster  than  the  IFDPE  model, 
especially  for  problems  involving  sloping  bathymetry. 
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